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Abstract 



When highly active antiretroviral therapy is administered for long periods of time to HIV-1 
infected patients, most patients achieve viral loads that are "undetectable" by standard 
assay (i.e., HIV-1 RNA < 50 copies/ml). Yet despite exhibiting sustained viral loads be- 
low the level of detection, a number of these patients experience unexplained episodes of 
transient viremia or viral "blips" . We propose here that transient activation of the immune 
system by opportunistic infection may explain these episodes of viremia. Indeed, immune 
activation by opportunistic infection may spur HIV replication, replenish viral reservoirs 
and contribute to accelerated disease progression. In order to investigate the effects of con- 
current infection on chronically infected HIV patients under treatment with highly active 
antiretroviral therapy (HAART), we extend a simple dynamic model of the effects of vacci- 
nation on HIV infection [Jones and Perelson, JAIDS 31:369-377, 2002] to include growing 
pathogens. We then propose a more realistic model for immune cell expansion in the pres- 
ence of pathogen, and include this in a set of competing models that allow low baseline 
viral loads in the presence of drug treatment. Programmed expansion of immune cells upon 
exposure to antigen is a feature not previously included in HIV models, and one that is 
especially important to consider when simulating an immune response to opportunistic in- 
fection. Using these models we show that viral blips with realistic duration and amplitude 
can be generated by concurrent infections in HAART treated patients. 



Introduction 



Adherence to a regime of highly active antiretroviral therapy (HAART) suppresses the viral 
loads of most chronically infected HIV patients below the level of detection by standard as- 
say. However, a number of these otherwise well-suppressed patients experience unexplained 
"viral blips" while on therapy. Di Mascio et al. (2003a) in a study of 123 patients found 
that the mean blip amplitude was 158 db 132 HIV RNA copies/ml, with the distribution 
skewed towards low amplitude blips. It addition, Di Mascio et al. (2003a) suggest that a 
viral blip is not an isolated event but rather is a transient, intermittant episode of detectable 
viremia (HIV-1 RNA > 50 copies/ml) with a duration of roughly two to three weeks They 
further showed that the amplitude and frequency distribution of these viral transients are 
consistent with viral load rising sharply, followed by slower, two-phase (double exponential) 
decay (Di Mascio et al., 2003a). An example of a typical blip is shown in Figure ^ A. The 
frequency of these episodes appears to be inversely correlated with CD4+ T-cell count at 
baseline, prior to initiation of drug therapy (Di Mascio et al., 2003a, 2004b). Figure 
shows data from the 123 patients in the Di Mascio et al. 2003a study. CD4+ T cell counts 
at the onset of therapy are plotted against the frequency of blip occurances per sample, 
showing a positive correlation between low initial CD4+ counts and blip frequency. 
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Figure 1: Generalized viral blip shape, and correlation with CD4+ T-cell count at the onset of 
therapy. A, stylized blip with rapid rise and two-phase exponential decay. Given a mean duration 
over threshold of 20 to 30 days, the two decay constants assume values of roughly 0.6 and 0.06<i -1 
(Di Mascio et al. 2003a), consistent with the mortalities of short lived productively infected cells 
and longer-lived chronically infected cells. B, CD4+ count at the start of therapy versus viral blip 
frequency. Both panels after Figure 4 (panels B,C), Di Mascio et al. 2003a. 
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A number of possible causes of blips have been suggested, including but not limited to missed 
drug doses, activation of latently infected cells and consequent release of virus, release of 
virus from tissue reservoirs, and a rise in target cell availability due either to vaccination 
(Jones and Perelson, 2002) or one or more coinfections by opportunistic pathogens, which 
then increase viral replication. The observation of Di Mascio et al. (2003a) that blip 
frequency is inversely correlated with baseline CD4+ T cell count suggests that patient 
specific factors, such as susceptibility to infection, which increases at low CD4 counts, 
may play a role in blip generation. Prior work on untreated, chronically infected HIV 
patients documents increases in viral load associated with vaccination (O'Brien et al., 1995; 
Staprans et al., 1995; Brichacek et al., 1996; Stanley et al., 1996) and with opportunistic 
infection (e.g., Donovan et al., 1996). McLean and Nowak (1992) proposed some of the first 
models of enhanced HIV replication due to immune stimulation via opportunistic infection, 
and showed how the positive feedback between enhanced HIV replication and incomplete 
immune control of pathogens due to HIV-immunosuppression leads finally to immune failure 
and full-blown AIDS. 

This paper examines the hypothesis that viral blips result from random encounters with 
replicating antigens - or transient opportunistic infections in HIV positive patients. We be- 
gin with a simple model for coinfection, explore its biological shortcomings, and develop a 
series of incrementally more complex models. To minimize the introduction of new param- 
eters while accurately simulating the dynamics of transient viremia, we add missing biology 
at each step and then test each incremental model to see where it falls short and where it 
is sufficient. In the process, we introduce biologically realistic mechanisms for programmed 
immune cell proliferation, and include features which allow robust low viral loads under 
drug treatment. 

A simple model for coinfection 

We begin with the following simple model for coinfection in the presence of HIV infection, 
which is a generalization of the HIV infection model developed by Perelson et al. (1996, 
1997) and reviewed by Perelson (2002). 
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Here the antigen, A, is a growing pathogen, T are uninfected CD4+ T cells, T* are cells 
productively infected with HIV, C are cells chronically infected with HIV, and V represents 
HIV-1 (RNA copies/ml). Pathogen A undergoes density-dependent growth described by a 
logistic law with maximum growth rate ro and carrying capacity A max . As in earlier work 
by Jones and Perelson (2002), we assume that the antigen is cleared in a T cell-dependent 
manner with rate constant 7. We further assume that uninfected T cells, T, are generated 
at rate A, die at rate dr, and are infected by virus with rate constant k. Assuming reverse 
transcriptase (RT) inhibitors are administered, the infectiousness of the virus k is reduced 
by (1 — e), where e is the efficacy of the RT inhibitors and < e < 1. 

In this simple model we assume T cells are activated into proliferation at a maximum rate 
a in the presence of pathogen, and that the proliferation rate depends on the pathogen 
concentration with a half-saturation constant K. HIV infection of T cells results in pro- 
ductively infected cells T*, which die at a rate 5, and chronically infected cells C, with 
mortality fi < 5. A fraction of infection events a « 1 results in chronic infection. Chroni- 
cally infected cells live much longer, producing virus more slowly than productively infected 
cells. The inclusion of the chronically infected pool is motivated by the suggested two-phase 
decay of a viral transient (Di Mascio et al., 2003a). 

Finally, virus is produced by productively and chronically infected cells at rates NtS and 
Ncn, respectively, where Nt and N c are average burst sizes for productively and chronically 
infected cells. Virus is cleared at a constant rate c per virion. 

Based on previous work we take as a typical set of parameter values A = 1 x 10 4 ml -1 , 
k = 8 x 10~ 7 ml/d, a = 0.195, N T = 100 and N c = 4.11 (Callaway and Perelson, 2002); 
d = 0.01 d- 1 (Mohri et al., 1998); 5 = 0.70 d^ 1 and fi = 0.07 d^ 1 (Perelson et al., 1997); 
and c = 23 d -1 (Ramratnam et al., 1999). The antigen or pathogen clearance rate constant 
7 is a "fitted" parameter, which we set to 1 x 10 -3 ml d _1 . In a prior study of the effects on 
viral load of vaccination (with a non-replicating antigen) , this fitted value varied widely from 
patient to patient with values roughly 1 x 10~ 5 to 1 x 10~ 8 ml d _1 , reflecting differences in 
patient immune response (Jones and Perelson, 2002). We assumed a higher clearance rate in 
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this study since immune activation by a proliferating pathogen should be much greater than 
activation due to a small dose of non-replicating antigen. Note that where these parameters 
are used elsewhere in this paper, they retain the above values unless otherwise noted. 

In the absence of pathogen, and assuming chronic HIV infection, the viral load stabilizes at 
the following equilibrium state: 



V = -[{!- a)N T + aN c ] - ^— - , (2) 

implying an inverse relationship between steady state viral load and drug efficacy e. 

This model without the inclusion of drug therapy, i.e. e = 0, and with a non-growing 
pathogen, i.e., r$ = 0, worked surprisingly well for modeling HIV dynamics in chronically 
infected, untreated patients following vaccination with a common recall antigen (Jones and 
Perelson, 2002). 

However, this model does not generate realistic blips when a growing pathogen is substituted 
for a vaccine. There are several reasons for this: even in untreated patients, according to 
our model, T-cells respond very rapidly, eliminating the pathogen before it has time to 
grow (Figure ?? A), so there is only a relatively small immune response to the presence of 
the pathogen, and little change in viral load (Figure ?? B). The system remains relatively 
unaffected by increases in pathogen growth rate until ro reaches a critical point where 
the model immune system cannot completely eliminate the pathogen, and then there are 
predator-prey cycles. [Note that inclusion of a logistic growth term means that the pathogen 
can reach a carrying capacity, A max , but will not experience runaway growth.] With the 
addition of drug therapy (Figure ?? C) of efficacy e = 0.6461275 (see Appendix), baseline 
viral load is suppressed to 25 RNA copies (ml -1 ) and opportunistic infection results in a 
small, slow rise and fall in viral load, rather than a burst of viremia. If baseline viral load 
is suppressed further, then there is no appreciable change in viral load. 



Building a new model 



While different parameters could be explored, we believe that there is a need to generate 
more biologically realistic models. First, the majority of patients with viral loads below 50 
copies/ml who have been examined with more sensitive assays have viral loads that are still 
detectable, i.e. been 1 and 50 copies/ml (Dornadula et al.,1999, Di Mascio, 2003b). This 
sugests that appropriate models must exhibit robust low viral steady states for patients 
on therapy. By contrast, the model given by eqs. (^Q) requires drug efficacy e precisely 
tuned to many decimal places in order to yield a low, yet not vanishingly small, viral 
steady state (Callaway and Perelson, 2002). Second, recent immunological data suggest 
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Figure 2: Simulations of a simple coinfection model Q for pathogen growth rates tq — 1,2,3,4 
day" 1 (solid, dashed, dotted line in each panel) given no drug treatment (panels A, B) and treatment 
with a drug of efficacy e = 0.6451275 (panel C). A, Antigen growth assuming e = and Aq = 1; B 
Viral load (HIV-1 RNA/ml) assuming e = 0; C, Viral load (HIV-1 RNA/ml) assuming e = 0.6461275. 
Baseline VL is 25 RNA copies/ml. 



T-cell proliferation stimulated by a pathogen involves a cascade of divisions, all triggered 
by a brief exposure to antigen. Models should thus account for this type of "programmed" 
response. Lastly, immune response to a pathogen involves generation of effector cells that are 
responsible for clearing the pathogen. For many viral infections, CD8+ T cells are critical for 
containing and clearing the infection (Wong and Palmer, 2003). This is well documented 
in lymphocytic choriomeningitis virus (LCMV) infection, where a strong antigen-specific 
CD8+ T cell response is induced, leading to rapid elimination of the virus (Murali-Krishna 
et al., 1998a,b). 



T-cell proliferation under antigenic stimulation 

Upon exposure to antigen, naive CD8+ T cells undergo a burst of proliferation, entering 
a programmed cascade of divisions that culminate in the production of mature, activated 
effector cells (Kaech and Ahmed, 2001). This is followed by a programmed contraction in 
which most of the effector cells are subject to apoptosis, leaving a small, stable memory 
population (Badovinac et al, 2002). Revy et al. (2001) proposed a system of ordinary 
differential equations to analyze and describe T-cell proliferation. We adapt their model to 
incorporate activation by a growing pathogen, A, which is finally cleared by CD8+ effector 
cells, E. 
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Here again A is an opportunistic pathogen, tq is the pathogen growth rate, and 7 is the 
clearance rate constant for the pathogen, po and p are constant proliferation rates, Nq is 
the initial, naive cell pool, and the iVj are proliferative phases, or "division classes": for 
each i, the number of cells that have completed i divisions, and E are mature, pathogen- 
specific "effector" cells. As the CD8+ T cells proliferate they also differentiate into effectors, 
E. Here we assume cells become effectors after 4 divisions and stop proliferating after 8 
division; E = J2^=4^i- Alternatively, one could assume a fraction of cells, become 
effectors after i divisions, i.e., E = Yli=i Pi^i, or leave the maximum number of divisons 
as an adjustable parameter. Proliferative, non-effector phases (division classes A^o, N\, 
^3) undergo mortality at a rate d « d E , the death rate for activated effector cells. 

Experiments suggest that when quiescent cells are stimulated into proliferation, the initial 
cellular division takes longer than subsequent divisions (Gett and Hodgkin, 2000), and that 
the time to first division depends on features of the antigen stimulation. For example, using 
anti-CD3 antibodies rather than antigen to stimulate CD8+ T cells, Deenick et al. (2003) 
found that decreasing the anti-CD3 concentration lengthened the time to first division. We 
thus assume the rate of the first CD8+ T cell division depends on antigen according to a 
Type-Ill functional response, 



Po{a)=Po WTkJ) (4) 

where Kg represents the critical antigen level required to stimulate a CD8+ response [Table 
I]. Thus when A is low, i.e., A « Kg, growth is slow, whereas when A attains values greater 
than Kg rapid growth ensues, saturating at rate p$. Here n is a parameter frequently called 
the Hill coefficient and it determines the steepness of the response when A is near the 
"threshold" K 8 . 

While the initial division may be relatively slow, subsequent divisions occur more rapidly, 
without an initial delay, as cells are already activated. These divisions occur even if antigen 



8 



(A): T— cell stimulation 



(B): Division Classes 




Figure 3: Simulated results for CD8+ T-cell proliferation stimulated by a growing antigen. A, 
antigen [Aq = 1), shown in dashed red line; total CD8+ effector cells (solid black line). B, cell 
counts for the various division classes (dotted blue); total effectors (5^ i=4 2Vj) (solid black line), and 
naive cells (dashed green line). Since the process of proliferation is a cascade involving both growth 
and mortality and we assume growth is more rapid than mortality each successive division class is 
slightly larger than the preceding class. The effector population (black line) is subject to higher 
mortality than the less differentiated preceeding classes. 



is removed (Badovinac et al., 2002). We therefore assume subsequent divisions occur at a 
constant rate, p > po (Table 1). 

In the model we do not explicitly account for precursor cells that fail to divide, although in 
the experiments by Gett and Hodgkin (2002) such cells are present. 

Based on CD8+ T-cell response estimates for LCMV in mice (De Boer et al., 2001), we set 
p = 2.92 d _1 , and chose the initial activation/first division rate po = 1 d _1 . The initial 
quiescent cell pool is small, Nq = 500 cells/ml: it is estimated that one in 10 5 T cells is 
specific for any given epitope, and assuming that a healthy human being has roughly 10 11 
CD8+ T cells, we expect about 10 6 antigen-specific cells. Most naive cells are in lymphoid 
tissue, which in a 70 kg man weighs about 700 grams. Thus, with some cells circulating in 
the blood and others not naive, a density of about 500 cells/ml in lymphoid tissue, where 
most proliferation occurs, seems appropriate. 

Figure |21 shows results from simulation of equations (j3J). In this simulation, we assume eight 
division classes, Ni, i = 1, ...,8, plus naive cells, Nq, with effector cells, E, corresponding to 
the sum of the last four division classes. In panel A, the growing antigen (dashed red line), 
and CD8+ effectors (solid black line) are shown. Panel B shows the CD8+ T cells in the 
various division classes and the mature effector cells. Note that effector pool grows and is 
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depleted within the first 20 days, while the proliferating cells that have not differentiated 
into effectors are longer lived. 



Density-dependent infected cell death 

In the simple coinfection model Q, infected cells are assumed to have constant mortality: 
5 for productively infected cells, T*, and \i for chronically infected cells, C*, respectively. In 
this model under very effective drug therapy the steady state viral load becomes vanishingly 
small. However, if the death of infected cells is made density dependent, as one might expect 
when infected cell death is immune-system mediated and immune response is a function of 
the number of infected cells present, then low steady state viral loads are possible (Callaway 
and Perelson, 2002). 

Here we assume that the magnitude of the immune response against HIV-infected cells, and 
thus their death-rate, is an increasing function of the infected cell density (Callaway and 
Perelson, 2002). Holte et al. (2001), fitting data on viral decays in HIV treated children, 
suggested that this can be modeled by a power law, that is, by replacing the constant 5 
with a function 8*(T*) that depends on the total number of infected cells raised to a power 



s*(T; otal )=5'(T; otal r (5) 

and the constant \x with 

where u < 1 and T* otal = T* + C is the total number of HIV-infected cells, which represents 
the stimulus for a cell-mediated response to HIV- infection. 

In the following series of models employing density-dependent infected cell death, we explore 
various scenarios for CD8+ mediated immune responses to an opportunistic infection. 



Simplest case: An HIV infection model with programmed proliferation of 
CD8+ T cells 

We now combine the basic HIV infection and treatment model with CD8+ T cell response 
to pathogen in the form of the programmed proliferation model described above. Since 
the effector cells are pathogen specific, not HIV specific, the dynamics of the CD8+ T cell 
population remains decoupled from the HIV infection model, except through their indirect 
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effect on the pathogen. To the modef for programmed proliferation of effector cells under 
antigenic stimulation (equations I3J) we add the following model for HIV infection with 
density dependent mortality of infected cells: 
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Target cells, T, are activated by antigen following a type-II functional response, as in the 
basic model 



and are susceptible to infection by HIV at a rate k. We do not separate out antigen-specific 
CD4+ cells; instead we let the activation parameter a incorporate the fact that only a 
fraction of CD4+ cells in that pool may be activated by antigen. 

Because the model includes density-dependent cell death, we decouple viral production 
rates from cell death rates, and assume virus is produced at constant rates q\ = NtS = 70 
d _1 and q 2 = 0.28 d _1 , from short-lived and chronically- infected cells, respectively. Rate 
constants and parameters describing the proliferating CD8+ pool remain as above. A listing 
of additional fixed parameter values can be found in Table 1. 

Figure |I] shows that this model can generate transient viral load increases. Panel A shows 
that the pathogen initially grows and then is cleared (dashed line) as the effector cells (solid 
line) increase. Pathogen growth and clearance rates are variable, but for this example 
we assume r$ = 2 d _1 , and clearance rate constant 7 = 1 x 10 -3 . The effect of this 
pathogen growth is to stimulate increased infection and viral production (Fig. E£$). As drug 
treatment is assumed, CD4+ levels can be quite high, leading to an unrealistic number of 
activated target cells, and excessively large increases in viral load. Thus while the model 
does indeed produce viral transients, there is a tradeoff between the height and the width of 
the resulting "blip" . That is, one can obtain an appropriately narrow spike in viremia but 
it is unrealistically large, or one can obtain a blip of appropriate amplitude (i.e., 100-1000 
copies/ml, cf. Di Mascio et al. 2003a) but of unreasonable duration (upwards of 100 days). 

In our simulations, parameters associated with steady-state baseline viral load are fixed as 
discussed above and in the Appendix, and parameters associated with pathogen growth 
and CD4+ activation are considered variable. Results from sensitivity testing of the model 
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Figure 4: Simulations and sensitivity testing of equations I|tjl4l7|l . A, Antigen growth (dashed line) 
and CD8+ effector cells (solid line); B, Viral load (RNA copies/ml, in solid line), productively 
infected cells (ml' 1 , dashed line), and chronically infected cells (ml -1 , dotted line); C, Viral load 
(copies/ml) from sensitivity testing of CD4+ T cell activation parameter a, as a assumes the values 
a = 0.25,0.50,0.75, 1 (solid, dashed, dash-dotted and long-dashed line, respectively), ro = 2 in this 
example. D, Viral load (copies/ml) from sensitivity testing of pathogen growth K4, as K4 assumes 
the values K4 = 1, 10 1 , 10 2 , 10 3 d _1 (solid, dashed, dash-dotted and long-dashed line, respectively). 
Note that activation parameter a is held at 0.5 when not varied. 
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(varying of parameters associated with activation of target cells in response to pathogen 
growth and especially the growth of the pathogen itself) are shown in Figures |lp and D. 
The model is very sensitive to changes in the CD4+ T cell activation rate a: between the 
values a = 0.25 d~ x and a = 0.75 d _1 , the model goes from producing no viral transient to 
producing an unreasonably large viral transient (Figure0P). Variation in the half-saturation 
of pathogen load required for CD4+ stimulation, i<Q, shows the opposite trend: as the value 
of K4 and thus the pathogen load required for CD4+ stimulation increases, resulting viral 
transients go from being excessively large (Figure EJ), red line), to non-existent. 

Variation in tq shows that beyond a threshold value there is relatively small variation in 
blip peak amplitude with substantial, though biologically reasonable, variation in pathogen 
growth rate (results not shown). 

Interestingly, this model exhibits a delay between the peak in opportunistic infection and 
the peak in viremia (Figure Ef3). The delay appears to be related to the time it takes for 
the T cell pool to increase and generate a rise in the numbers of infected cells, T* and C, 
which then in turn generate more virus. As one sees (in Figure 0P), increasing the CD4+ 
T cell activation rate causes the peak in viremia to occur earlier. Whether these delays 
are realistic is unknown, but because antigen-specific CD4+ cells and "non-specific" CD4+ 
cells are all in the same pool and the pool has the relatively low mortality appropriate for 
"resting" cells, the T cell response may not be realistically modeled. 

Including separate antigen-specific and non-specific target pools 

In the prior model, the entire CD4+ pool was activated by antigen, resulting in an overly 
large target population and correspondingly large viral transients. To make the model 
more realistic, we now divide the target (CD4+) cells into non-specific and antigen-specific 
pools, T and Ta, respectively, where only the antigen-specific CD4+ T cells are activated 
by antigen. We couple the model for antigen-induced response of CD8 effectors with 
the following model for HIV infection of antigen-specific and non-specific CD4 cells: 
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uX + f(A)T A - (1 - e)k A VT A - d TA T A (8a) 

(l-e)k A VT A -(6*(u))TX (8b) 

(8c) 

(1 - u)X - (1 - e)WT - d T T (8d) 
(l-e)(l-a)jfeVT-((5*(w))T* (8e) 
a{l-e)kVT-n*{u)C (8f) 
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where £ = £]f =4 iVj, and p (^) is defined in (@J) and /(^4) defined in Q). 

Antigen-specific CD4+ T cells, Xa, are assumed to be small fraction v of the CD4+ pool T 
(see Table 1), and are activated by antigen following a type-II functional response, as in 0. 
The value of v is poorly defined, and may vary widely from roughly 10~ 5 to 10~ 2 (percent 
activated 0.001% to 0.1%, respectively). Response to a particular epitope generally occurs 
at a frequency of about 10~ 5 (0.001%), but a pathogen has many epitopes, thus as many as 
10 to 20 responses might occur for a single pathogen, yielding a frequency of 10~ 4 (0.01%). 
However, if one considers levels of activation in HIV-infected versus healthy people, one 
finds that in HIV-infected people, roughly 6.5% of CD4+ T cells are activated, depending 
on total CD4+ count, versus about 1% in healthy individuals. This level of activation 
might reflect responses to a number of pathogens, or in the case of immuno-compromised 
HIV patients, may be a homeostatic response to low CD4+ T cell counts. Based on these 
considerations, we choose a frequency v = 0.01, corresponding to an activated fraction of 
about 1%. (Note that we could easily choose a slightly higher or lower fraction and obtain 
similar results by adjusting the activation parameter, a 0). 

Once activated, target cells are subject to infection by HIV at an increased rate k A . The 
subsequent increase in viral load results in increased infection of the non-specific pool, 
T. Both infected antigen-specific T* A and infected non-specific cells T* are cleared in a 
density dependent fashion. Rate constants and parameters for the proliferating CD8+ pool 
remain as above. Analysis of SIV infection of resting and infected CD4+ T-cells suggests 
that activated, infected CD4+ cells produce on average twelve times more virus than non- 
activated infected cells (Zhang et al., 2004). We thus set q% = 12gi in equations 0- 

Using density-dependent infected cell death, the model can produce realistically low baseline 
viral loads under treatment. For this example the baseline viral load is about 3 RNA copies 
per ml. Opportunistic infection results in detectable viral transients (Figure ISJA.) , with the 
stimulus for induced replication arising from the (small) antigen-specific CD4+ T cell pool, 
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A: Viral Load 



B: Infected CD4+ target cells 




Figure 5: Simulation of equations JHJ for a patient assumed to be in steady state with V = 3 
copies/ml. At t = a pathogen (initial value Aq — 1), growing at rate tq = 2.0 days" 1 , is 
encountered and eliminated by an immune response with rate constant 7 = 1 x 10~ 3 . Here v = 0.01, 
e = 0.62, a = 1.2, Kg = 1000, K4 = 100. A, Viral load (copies/ml). Dashed line shows limit of 
assay detection at 50 (copies/ml). B, Non antigen-specific productively infected (solid line) and 
chronically infected (dashed-dotted line) CD4+ cells, as well as antigen-specific infected CD4+ cells 
(dashed line), increased by a factor of 10 to better show dynamics. 
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but most of the increase in viral load coming from infection of antigen non-specific T cells 
(Figure |S|3). The resulting transients have a sharp rise-time and a slow decay, consistent 
with observation (Di Mascio et al., 2003a). However, blip duration is generally too long, 
upwards of 30 to 90 days, respectively, for transients of reasonable amplitude (100-1000 
RNA copies per ml). 

In addition, the model is very sensitive to changes in the parameters a and K%, with less 
extreme variability in viral load resulting from variation of the pathogen growth parameter 
tq and the half-saturation parameter K4 (Figure EJ). Small variation in the CD4+ cell 
activation parameter a alone results in a spectrum of behavior from no detectable viremia 
to excessively large bursts in viral replication (Figure El A) . In addition, though in general 
infectivity k is poorly constrained, infectivity for activated antigen specific cells Ica must be 
set at levels fifteen times higher than that of quiescent cells to produce detectable viremia, 
and this may not be realistic. 

Effects of parameter variation 

• Varying the CD4+ activation parameter, a, results in large changes in peak viremia 
over a relatively small range of parameter values, as in Figure IHJA Note that blip 
durations (i.e., decay times) are too long in all of these examples, and only begin to 
become shorter with increased values of a (i.e., greater proportions of the transient 
arising from infection of activated, antigen-specific cells). 

• Increasing ro results in more aggressive pathogen growth and swifter CD8+ immune 
response. For lower values of ro, the CD8+ response is delayed, and f(A) (J7| remains 
elevated longer. A smaller ro thus results in a larger amplitude viral transient, with 
ro = 1 days -1 resulting in the largest blip, and successively larger values of ro resulting 
in smaller blip sizes (Figure ©3). 

• Increasing the half-saturation constant for CD8+ proliferation, K$ in has a positive 
effect on the sizes of the resulting viral transients, as this allows greater pathogen 
growth before an immune response begins, while increasing the half-saturation for 
pathogen growth in the activation function f(A) has a negative effect on blip size, 
as the slope of the "activation curve" for target cells becomes gentler (Figures |HP,D). 

A complex but more biologically realistic version of this model would allow proliferation of 
the antigen-specific target pool - yielding more activated target cells at the effector stage 
- and might also allow infection of the proliferating division classes. This would boost the 
amplitude of the transient and lower the value of k^ required to produce detectable viremia, 
thus addressing two flaws in the present model |HJ). We explore both of these options below. 
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Figure 6: Sensitivity testing, equations ©■ Unvaried parameters are as in Figure[S] A, Viral load 
as the activation parameter a varies from 0.5 (solid line) to 1.5 (dash-dotted line); B, Viral load as 
pathogen growth rate rg varies from 1 (solid line) to 4 (dash-dotted line); C, Viral load as Kg, (@J 
varies from 10 2 (solid line) to 10 5 (dash-dotted line); D, Viral load as K4 varies from 1 (solid 
line) to 10 3 (dash-dotted line). 
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Programmed proliferation of both CD8+ (effector) and antigen-specific 
CD4+ (target) cells 

CD4+ T cells as well as CD8+ T cells undergo programmed proliferation upon exposure to 
antigen (Lee et al., 2002). These proliferating cells are more susceptible to HIV infection 
than resting cells. Thus incorporating CD4 cell proliferation cascades should affect HIV 
replication. We incorporate CD4+ proliferation of antigen-specific CD4+ T cells, with each 
proliferating division class after the first division assumed activated and thus infectable by 
HIV-1. The model is identical to equations (|8|). with the addition of sets of equations for 
the proliferating To, Tj,... classes, and the infected T*, T 2 * classes. 

As in prior models we assume CD8+ cells undergo eight divisions, so there are eight division 
classes with the assumption of effector status after the fourth division. CD4+ T cells have 
slower kinetics than CD8+ T cells, and undergo fewer programmed cycles of proliferation. 
In the case of LCMV infection of mice this results in 20-fold lower expansion and a delay by 
one day until the peak of expansion is reached (De Boer et al., 2003). To limit introduction 
of new parameters, and because little is known about human CD4+ proliferation rates, we 
assume the same proliferation rates for CD8+ and CD4+ cells. Since CD4+ cells undergo 
fewer divisions, we assume four CD4+ division classes as opposed to eight, and further 
assume that they are activated by pathogen and infectable after the first division, and 
that they suffer elevated mortality 8a when activated. In addition, we allow continued 
proliferation (division) of infected division classes, which divide to produce infected cells. 
The initial number of naive cells in both the CD8+ (i.e., iVo) and CD4+ (i.e., To) pools is 
set at 500 cells. To the model for programmed proliferation of CD8 effectors under antigenic 
stimulation we add the following equations for proliferation and infection of antigen-specific 
CD4 cells and the passive infection of the nonspecific CD4 pool: 
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dT 

Antigen-specific CD4 = -(po(A) + d )T (9a) 

dT\ 

dTi hffi-i - (p + d)Ti - (1 - e)/c A ^,, i = 2, 3 

\2pT; - (1 - e)k A VT i+1 - d TA T i+1 ,i = 4 



2p (A)T - (p + d)Ti - (1 - eJfc^Ti (9b) 



(l-e)fc A Fri-(**H+p)3t (9d) 

(9e) 



cG? Ul-e)k A VT i + 2pT*_ 1 -(6*(Lo)+p)T*,i = 2,3 

1(1 - e^FZ} + 2pr f *_! - (<5*(w))2?,t = 4 



Non-specific CD4 — = A - (1 - e)£;FT - d T T (9f) 



dt 
dT* 

~df 
d£_ 

~dt 
dV_ 

~dt 



(l-a)(l-e)kVT-(5*(u))T* (9g) 
a(l-e)A;VT-/i*(w)C7 (9h) 
q x T* + q 2 C + qsJ^T? ~ cV (9i) 



where E 1 = Y1^=a an( ^ Po(A) is defined in equation (j2J. 

Figures [7] and show simulation results and sensitivity testing for equations © . Figure 
shows viral load, panel B shows productively and chronically infected cells (solid and 
dashed line, respectively), and panel C shows infected antigen-specific CD4+ effectors (bold 
solid line) and division classes (thin dashed line) for a representative parameter set (ro = 2 
d _1 , Kq = 10 3 , and k A = 6 x 10~ 6 ). The rise in infected and chronically infected cells 
from the non-specific pool is fairly small, but most of the new replication in the transient 
is coming from infection of the even smaller activated antigen-specific pool. The timing of 
the peak in new productively and chronically infected cells from the non-specific pool is too 
late to account for the bulk of the viremia in the transient. 

The model is fairly robust, producing viral transients of biologically reasonable amplitude 
and duration over a range of parameter values. In Figure |H1 we test the model for reasonable 
variation in three parameters critical to pathogen growth and viral replication. Results for 
all three tests are shown in linear scale, as the variation in peak viremia spanned only one 
order of magnitude. Panel A shows variation in viral load (copies/ml) as rg is varied from 
ro = 1 d _1 to ro = 4 d _1 , while parameters Kg and k A are held at the values listed above, 
and all other parameters are held constant as listed in Table 1. Unlike the prior model 
(equations |S1 and Figures [5] and |f)J) , in this case peak viral load is an increasing function 
of ro, though the increases in VL with each increase in ro are moderate. In panel B, ro 
is again set to ro = 2 d _1 , and the parameter Kg is varied from 10 2 (solid line) to 10 5 
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Figure 7: Simulation of equations (J5J. Again a patient is assumed to be at steady state (V = 3 
copies/ml) when at t = a growing pathogen is encountered. Here, fc^ = 4 x 10 _6 ml/d, r = 2.0 d _1 
and 7 = 1 x 10~ 3 ml/d. Other parameters are as described in text and shown in Table 1. A, Viral 
load (copies/ml); dashed line shows limit of assay detection at 50 (copies/ml). B, Total productively 
(solid line) and chronically infected CD4+ cells (dashed line), C, Infected antigen-specific CD4+ 
effectors, solid line, plotted with infected antigen specific CD4+ division classes, thin dashed line. 
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Figure 8: Sensitivity testing for equations Each panel shows viral load (copies/ml) on a linear 
scale. A, Viral load as pathogen growth parameter varies, ro = 1, 2, 3, 4 d _1 in solid, dashed, dotted, 
and dot-dash line, respectively; B, Viral load as K% © varies from I0 2 (solid line) to I0 5 (dash- 
dotted line); C, Viral load due to variation in antigen-specific, activated cell infectivity fc^ = 8 x 10~ 7 
(cf. Callaway and Perelson, 2002), 1.6 x 10~ 6 , 3.2 x 10~ 6 , 6.4 x 10~ 6 , 8.0 x 10" 6 (ml-KNA^ 1 ^ 1 ); 
solid, dashed, dotted, and dot-dash line, respectively. 
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(dash-dotted line), resulting once again in moderate increases in peak viral load and blip 
duration. Variations in the infectivity of the antigen-specific cells Jza over a biologically 
reasonable range of values also results in moderate and acceptable increases in peak blip 
viremia (Figure |HP). 



Models with heterogeneous target pools 

Another intuitive explanation for the presence of intermittent viremia is to assume that 
different cell types, or cell populations, might have different drug penetrances, perhaps due 
to different levels of expression of P-glycoproteins that have the potential to pump drugs 
out of cells (Kim et al., 1998). Having cell populations with reduced drug penetration may 
produce episodes of transient viremia in the presence of opportunistic infection. Indeed, 
there is evidence from in vitro studies of heterogeneity in intracellular drug concentrations 
(c/. Kim et al, 1998; Perno et al, 1998; Puddu et al, 1999). In the following model 
we include two co-circulating populations of target cells with differing drug penetration, 
wherein drug efficacy is reduced in one cellular population. This is, in addition, a simple 
means of producing low steady-state viral loads in the presence of drug therapy which does 
not invoke density dependent mortality of infected cells (c/. Callaway and Perelson, 2002). 

The model again includes expansion of antigen-specific CD8+ cells, with the effector class 
E as the final division class. Effector cells have the same mortality as mature CD4+ T-cell 
classes, but remain decoupled from the two co-circulating CD4+ classes T\ and T%. Target 
cells do not proliferate from the naive state but are assumed to have passive sources Ai 
and A2 (Callaway and Perelson, 2002). The model assumes differential drug efficacy in 
the co-circulating target cell populations: in one population (Ti), drug efficacy is e, and 
in the second population (T2), drug efficacy is fe, where efficacy has been reduced by a 
factor / < 1. Finally, the model includes two separate chronically infected cell pools, C*, 
derived from infected non antigen-specific T-cell pools, T* and Again, the model for 
programmed proliferation of antigen-specific CD8 effector cells is as in (J3J). 
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Figure 9: Simulation of equations l|10|) . Here, k\ = 8 x 10~ 7 , = 1 X 10~ 5 and other parameters 
are as described in text and shown in Table 1. A, viral load (solid line, in RNA copies/ml) assuming 
ro = 2.0 and 7 = 1 x 1CP 3 . B, Productively and chronically infected CD4+ cells, dashed and solid 
lines respectively. Note that these cells are contributed only by the cellular pool for which drug 
efficacy is assumed to be imperfect. C, Viral load (log±o) as the activation parameter a varies from 
0.01 (solid line) to 1 (dash-dotted line); D, Viral load (linear scale)as pathogen growth rate rg varies 
from Id -1 (solid line) to 4 (dash-dotted line); E, Viral load as K 8 (QJ varies from 10 2 (solid line) to 
10 5 (dash-dotted line); F, Viral load as K 4 J7J varies from 10 1 (solid line) to 10 4 (dash-dotted line). 



22 



-± = X 1 + f(A)T 1 -(l-e)k A VT-5T (10a) 
= X 2 + f(A)T 2 -(l-fe)k 2 VT-8T (10b) 

dT* 

-± = (1 - e )(l - a)kVTt - ST^ (10c) 
dT* 

= {l-fe){l-a)WT 2 -5TZ (lOd) 

dC 

-j± = a{l - e)WT 1 - nd (We) 
dC 

—1 = a(\-fe)WT 2 -iiC 2 (lOf) 
at 

dV 

— = N T S{T^ + 1^)+N c ii{C 1 + C2)-cV (lOg) 



where again E = 5Zj=4-^i> Po(A) is as in equation (J3J, and f(A) is as in equation (J7J). CD8 
effector cells are activated by antigen following a type-II functional response, as in the basic 
model ©• As in equations CD4+ cells in both pools are activated following a type-II 
functional response, with activation parameter a. Other parameters are as defined above 
and in Table 1, with the exception of parameters specific to a given target cell pool, which 
are given subscripts appropriate to that pool. Infection rates k A and k 2 and cell production 
rates Ai and A2 are thus specific to target cell pools T\ and T 2 , respectively. As in Callaway 
and Perelson (2002), we do not introduce different infected cell death rate constants 5 and 
[i for the differing target cell types. Assuming that e = 1 in population 1, the steady state 
viral load is 



V(e = 1) = %1 - a)N T + aN c ] - a \ . (11) 

Further assuming that V = 100 when / = 0, we set A2 = 57 cells (ml -1 d -1 , see Appendix) 
and k 2 = 10~ 4 ml copies -1 d _ , again as in Callaway and Perelson (2002). 

In Figure are shown simulations and sensitivity testing of ()10|) . Panel A shows viral 
dynamics assuming a pathogen with growth rate ro = 2; viral load is shown in solid line 
and the threshold for detection [50 copies/ml] is shown dashed. In Figure EJlB is shown the 
rise in productively and chronically infected cells for both target cell populations. Infected 
cells for target cell population 1 do not rise at all, so all new viral replication and thus 
the entire viral transient comes from newly infected cells in population 2. Finally, the 
above example shows that opportunistic infection will produce viral transients of typical 
size and duration in numerous models with suitable mechanisms for generating robust low 
steady-state viral loads. 
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Discussion 



Models of HIV infection have generally not included detailed models of immune responses. 
When immune responses are included in such models, they tend to be very stylized, phe- 
nomenological rather than mechanistic, and with simple predator-prey type dynamics (cf., 
Nowak and Bangham, 1996; Wodarz et al., 2000; Wodarz and Nowak, 2000). Recent exper- 
imental work has shown that both the CD4+ and CD8+ classes of T cells undergo bursts 
of proliferation on encounter with antigen, and that programmed proliferative responses 
continue even if antigen is removed (Badovinac et al., 2002). This is a feature that has 
not previously been included in HIV models, and one that is especially important to con- 
sider when simulating the effects of coinfection - infection with a growing pathogen - on 
chronically infected HIV patients. 

Here we have built a set of models of increasing realism that captures the proliferative 
response of both CD8+ and CD4+ T cells, and the dynamics of chronic HIV infection 
in the presence of opportunistic infection. In order to tie our models to data we have 
focused on explaining the rare occurrences of transient episodes of viremia, or viral blips, 
that are seen in patients treated with potent antiretroviral therapy. Blip frequency is 
independent of time on treatment, and if viral loads were to continue to fall towards 
extinction as they do in most standard models, one would expect the occurrence of 
detectable viral transients to become rarer with increased time on treatment, something 
which is not observed (Di Mascio et al., 2003a). Furthermore, chronically infected patients 
are usually on therapy and their viral loads are generally suppressed below, or close 
to, the level of detection by standard assay. For these reasons we focused our analysis 
on HIV models that allow robust low steady state viral loads (Callaway and Perelson, 2002). 

Simulation results 

Starting from a very simple model of T cell and HIV dynamics under drug therapy (equa- 
tions ^) , which when challenged with a growing pathogen only produced a gradual increase 
in viral load rather than a blip (Figure 1C) we built a family of new models, at each step 
adding more biological realism and complexity. The results from each successive interim 
model guided us towards the creation of biologically realistic and robust end-models, which 
are capable of producing viral transients of "typical" duration and amplitude, under a 
biologically reasonable parameter regime. 

Addition of a proliferating effector pool and density dependent infected cell death to the 
initial model results in a model which can produce viral transients, but with an inherent 
trade-off between blip duration and amplitude (equations EJ) • The model either produces 
a blip of appropriate duration, but of amplitude several orders of magnitude larger than 
desired, or of appropriate amplitude, and upwards of 100 days duration (Figure^} depending 
on the choice of parameters. This appears to be due to the fact that in the model the entire 
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pool of CD4+ cells is activated by antigen, resulting in a large spike in viral replication. 
In addition, this large target pool has the low mortality assumed for "resting" cells, which 
affects the duration of the resulting blips. 

We then separate the target pool into antigen specific and non-specific CD4+ T cells, 
and allow activation by antigen - and thus increased susceptibility to infection - of the 
antigen-specific pool (equations [SJ). Simulations of this model can produce blips with a 
sharp rise-time and slow decay, consistent with observation, and of "typical" amplitude as 
described above. However, even the smaller (< 100 RNA copies/ml) blips are slightly too 
long in duration ( 30 days) (Figure EK), and the model is very sensitive to changes in the 
activation parameter a and the half-saturation constant K$, both of which are difficult to 
quantify via experiment or observation. 

Allowing the antigen-specific target cell pool to proliferate upon exposure to antigen is a 
natural next step, which eliminates the biologically ambiguous activation parameter, a, and 
produces a model which generates "blips" of reasonable amplitude and duration over a 
range of parameter values (equations E0 Figure 

A model that includes cell populations with different drug susceptibilities (equations HOj) 
generates robust low steady state viral loads without having to invoke density-dependent 
cell death and its attendant parameters (cf. Callaway and Perelson, 2002). This model 
when coupled with a CD8+ T cell model of the response to a replicating antigen, can also 
produce viral transients of appropriate amplitude and duration. In this case it is interesting 
to note that the viral transient arises entirely from the second compartment with reduced 
drug penetration; thus the size and duration of such transients might scale with the size of 
the compartment from which they arise. 

Limitations 

The models presented here, despite their increasing complexity, remain simplifications of 
an in vivo system. In order to limit the size of our models and the number of parameters 
therein, we have not included HIV-specific CD8+ and CD4+ T cells. In addition, we have 
focused entirely on a T-cell mediated immunity, and have not included antibody (B-cell) 
responses. 

The clonal expansion and survival of CD8+ T cells may depend on the activity and "help" 
of CD4+ T cells, but mostly during formation of memory and during antigen rechallenge 
(Bevan, 2004). Thus primary CD8+ T cell responses are undiminished in a CD4+ T cell 
deficient system: when antigen alone invokes an overwhelming inflammatory response, as in 
many infections, primary CD8+ T cell responses are entirely independent of the CD4+ T 
cell response (Bevan, 2004). Yet CD8+ T cell memory during antigen rechallenge does not 
function properly in a CD4+ T cell deficient system. Furthermore, there is some evidence 
that frequency of viral transients is statistically correlated with low CD4+ counts (Di Mascio 
et al., 2003a). It would be worth exploring whether, in a CD4+ T cell deficient environment, 
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insufficient CD8+ T cell responses to chronic immune rechallenge contribute to episodes of 
transient viremia. While we do not include CD4+ "help", or indeed, CD8+ T cell memory 
in our system, these features would be necessary for any biologically realistic simulation 
of an immune system encountering pathogens in a stochastic fashion. We have omitted 
these features here as we have not modeled antigen rechallenge and immune memory. Our 
future research interests include stochastic exploration of immune activation during random 
encounters (and re-encounters) with pathogen, and thus will invoke immune memory. 

Conclusions 

We have shown that episodes of transient viremia in otherwise well-suppressed, chroni- 
cally infected HIV patients under drug therapy may be triggered by randomly- occurring 
opportunistic infections, which cause a rise in activated target cells, and thus transient 
bursts of infection and resultant viremia. Simple models for coinfection, which incorporate 
stylized, predator-prey type immune dynamics, cannot account for the occurrence of viral 
transients. However, if biologically realistic T cell dynamics that include programmed pro- 
liferation and contraction of the CD8+ T cell are incorporated into the model, along with a 
mechanism for producing robust low steady-state viral loads under HAART, then we have 
shown that opportunistic infection can serve as a "forcing function" for transient episodes of 
viral replication and can explain the observation of viral blips in otherwise well-suppressed 
HIV patients. This is just one plausible mechanism among many for the generation of blips; 
the list of possible causes also includes missed drug doses, activation of latently infected 
cells and release of virus from tissue reservoirs. It remains to be shown that these other 
mechanisms also work in a model and produce transients of appropriate amplitude and 
duration. 



Acknowledgements 

This research was supported by the Department of Energy under contract W-7405-ENG-36 
and by NIH grants AI28433 and RR06555 (ASP). We thank Michele Di Mascio for the data 
used in Figure 1. We thank Rob DeBoer for many useful discussions on the topic of T-cell 
proliferation, and the Cornell EEB eco-theory lunch-bunch, in particular Stephen Elmer, 
for comments which greatly improved this paper. 



26 



Appendix 



All simulations were written in the R programming language (version 1.9) and run on a 
Windows XP platform, using the odesolve package. 

Steady state values were obtained numerically by setting fixed parameter values and running 
the model until the simulation reached steady state. Steady state values for each variable 
were noted and used as initial conditions for simulations of the effect of concurrent infection. 

Parameter values: 

e: for the chronically infected cell model we seek a steady state viral load V of around 
50 copies/ml. Solving the expression for steady state viral load (J2J) for e, we obtain: 

dc/n 

e ~ ~ X[(l-a)N T + aN c ] -cV" 

Substituting in the appropriate parameters, one obtains a required drug efficacy of about 
0.6458 for a steady state viral load of V = 50. This efficacy is close to the critical efficacy 
for the models presented here, given our parameters. 

A2: In the two-compartment co-circulating cell model (|10|) . we invoke the expression for 
steady-state viral load, assuming e = 1 in the first cellular pool Q11JI - Solving for the 
parameter A2: 

, = VK 2 (l-f) + d2 

2 C K 2 (l-f)[(l-a)N T + aN c } 

Further assuming that V = 100 when / = and substituting in appropriate parameter 
values gives a value A2 = 57. 

fj,': Observing that the ratio <5 : fx = 10 for the death rates of productively versus chronically 
infected cells (Callaway and Perelson, 2002) we retain the same ratio of death rates under 
the assumption of density dependent mortality for both cellular pools. Thus, if 5' = 0.7863 
cr^mZ-ceir 1 )", then // = 0.07863 cT^mZ-ceir 1 )". 
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Table 1. Model Parameters 



Parameter 


Description 


Value 


Reference 


a 


T-cell activation parameter 


variable 




a 


fraction chronically infected 


0.195 


CP02 


Amax 


pathogen carrying capacity 


10 s (units/ml) 




C 


clearance rate 


23 (O 


R99 


d 


death rate, quiescent cells 


0.01 (O 


DB01 


d 


death rate, division classes 


0.1 (O 


Set 


d TA 


death rate, antigen specific 


0.325 (O 


DB03 


dx 


death rate, non-specific 


0.01 (O 


DB01 


5 


death rate, infected cells 


0.7 (d- 1 ) 


P97 


5' 


density dependent mortality 


0.7863 


CP02 






d- 1 {mbcc\\- 1 ) 0J 




d E 


death rate, effector cells 


0.325 (d- 1 ) 


DB03 


e 


drug efficacy 


< e < 1 


Appendix 


7 


pathogen clearance rate 


10- 3 (O 


see text 


k,ki 


infectivity, single(first) target pool 


8 x 10~ 7 


CP02 






(mZ-RNA- 1 ^ 1 ) 




k A 


infectivity, Ag-specific target pool 


6 x 10~ 6 


set 


k 2 


infectivity, second target pool 


10~ 4 


CP02 


K,K 8 


antigen half-saturation for stimulat- 


variable 






ing CD8 cells 






Ki 


antigen half-saturation for stimulat- 


variable 






ing CD4 cells 






A,Ai 


passive T cell source 


10 4 (cells /ml/d) 


CP02 


A 2 


passive T-cell source 


56 (cells/mZ/d) 


Appendix 


A*c 


mortality, chronically infected 


0.07 


P97 


M' 


density-dependent mortality 


0.07863 


Appendix 


N Q ,T Q 


initial quiescent population 


500 cells 


see text 


N c 


burst size, chronically infected 


4.11 RNA copies 


CP02 


N T 


burst size, productively infected 


100 RNA copies 


CP02 


V 


antigen-specific fraction 


0.1 


S98 


Pt) 


initial proliferation rate 


1 (d- 1 ) 


DB01 


P 


proliferation rate, classes 1, .., k 


2.92 (d- 1 ) 


DB01 


OJ 


See equation (0) 


0.01 


CP02 


91,92,93 


HIV production, density dependent 


70,0.28 (d- 1 ) 


CP02 


93 


HIV production, density dependent 


840 (d' 1 ) 


Z04 



CP02: Callaway and Perelson 2002, DB01: De Boer et al. 2001, DB03: De Boer et al. 
2003, P97: Perelson et al. 1997, S98: Sachsenberg et al. 1998, Z04: Zhang et al. 2004. 
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